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ABSTRACT 

A two-step hybrid perturbation-Galerkin technique for improving the usefulness of per- 
turbation solutions to partial differential equations which contain a parameter is presented 
and discussed. In the first step of the method, the leading terms in the asymptotic expan- 
sions) of the solution about one or more values of the perturbation parameter are obtained 
using standard perturbation methods. In the second step, the perturbation functions ob- 
tained in the first step are used as trial functions in a Bubnov-Galerkm approximation. This 
semi-analytical, semi- numerical hybrid technique appears to overcome some of the draw- 
backs of the perturbation and Galerkin methods when they are applied by themselves, while 
combining some of the good features of each. The technique is illustrated first by a simple 
example. It is then applied to the problem of determining the flow of a slightly compressible 
fluid past a circular cylinder and to the problem of determining the shape of a free surface 
due to a sink above the surface. Solutions obtained by the hybrid method are compared with 
other approximate solutions, and its possible application to certain problems associated with 
domain decomposition is discussed. 
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1. Introduction 


A two-step hybrid analysis technique, which combines perturbation techniques with the 
Galerkin method, has been presented and discussed by the authors. It was applied to some 
singular perturbation problems in slender body theory [5], as well as to several classes of 
problems involving ordinary differential equations [2,6,7]. That technique also promises to 
be useful in the analysis of a very wide variety of partial differential equation type problems 
as well. In this paper we apply the method to some problems involving several independent 
variables. In particular, we demonstrate its usefulness by applying it to a linear boundary 
value problem, a nonlinear boundary value problem, and to a (nonlinear) free boundary 

value problem. 

The method is based upon a hybrid technique which was apparently first studied by 
Ahmed K. Noor and collaborators in conjunction with the finite element analysis of geomet- 
rically nonlinear problems in structural mechanics (see Geer and Andersen [5] for several 
references). The Galerkin method has, of course, been known and used for a long time. 
However, a principle problem associated with its successful application lies in the choice of 
appropriate basis functions. In a series of papers Noor and his collaborators have shown for 
a variety of structural mechanics problems that the first few terms in a Taylor series expan- 
sion of the solution of a parameterized system of discretized equations can be particularly 
effective as Galerkin trial functions (or basis vectors). Subsequently, the present authors 
[2,5-7] have shown that the terms in either a regular or a singular perturbation expansion 
of the solution are also effective trial functions. In particular, we have demonstrated that 
the “reduced- basis” solutions can be useful for significantly larger values of the expansion 
parameter than the Taylor series or singular perturbation solutions on which they are based. 
A treatment of the reduced basis method from a mathematical point of view is given by Fink 

and Rheinbolt [3]. 

Some general observations about the technique are the following: 1) In many perturbation 
problems, much effort has to be expended to compute (analytically) each additional term 
in a perturbation expansion. Through the use of the proposed hybrid method, the known 
perturbation terms can be exploited more fully. 2) Another way of viewing the technique is to 
recognize that in many perturbation expansions the functional form of the higher-order terms 
can be well approximated by a linear combination of the lower-order terms. Thus, much of 
the effect of the higher order terms may be included by applying the reduced basis technique 
to a small number of lower order terms. 3) Since it is often possible to develop formal 
perturbation expansions for the solution about two or more values of the parameter, e.g. for 
small or large values of the parameter, the proposed method is a convenient way of combining 
the information contained in these different expansions. This is accomplished by allowing the 
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set of basis functions to include terms from the different perturbation expansions, which may 
be either regular or singular expansions (see [7] for several examples). 4) While the use of a 
Taylor series expansion is frequently limited by a finite radius of convergence, the proposed 
hybrid method can sometimes yield good results even well outside the radius of convergence 
(see [2] as well as our first example here), or even past a real singularity in the parameter [6j. 
In fact, when information from two or more expansions is employed, the method appears 
to provide meaningful (and often very accurate) approximations in “intermediate” regions 
of parameter values, e.g., in regions where the parameter is neither “small” nor “large” (see 
especially Geer and Andersen [7]). 5) The same general technique may also be applied in 
a fully numerical sense in that the perturbation functions themselves may be computed in 
discretized form. 6) However the perturbation functions are computed, the hybrid solutions 
share with the perturbation solutions the property that the accuracy is greatest near the 
expansion points and diminishes as the parameter moves away from these points. Usually 
a good indicator of accuracy is found by examining the difference in the hybrid solutions 
based onN-1 terms and N terms. Thus, it is possible to determine the range of parameter 
values for which the iV-term expansion is valid. 

In the following section, we describe the method in more detail and then apply it to a 
simple linear boundary value problem in Section 3. In Section 4, we apply the method to the 
problem of determining the flow of a slightly compressible fluid past a circular cylinder, while 
in Section 5 we apply it to the problem of determining the shape of a free surface induced 
by the presence of a sink above the surface. Each of the problems in sections 3-5 involves 
an elliptic boundary value problem, yet each serves as a model problem for a much wider 
class of problems. The problem of Section 3 is linear, and an arbitrary number of terms in 
its perturbation solution can be obtained in a straightforward manner. Thus, while an exact 
analytical solution is not known, it is possible to perform a rather detailed analysis of both the 
perturbation and hybrid results, including some explicit expressions for some of the hybrid 
solutions. In this problem the radius of convergence of the perturbation solution is finite 
and is limited by singularities which lie on the imaginary axis of the complex parameter 
plane. The problem of Section 4 involves a nonlinear partial differential equation, whose 
perturbation solution is tedious to obtain. Its radius of convergence is finite and appears to 
be limited by singularities which have real physical significance. The perturbation solution of 
the problem of Section 5 is the most difficult of the three problems to compute and is only an 
asymptotic (i.e., a nonconvergent) expansion [13]. For each of these problems we demonstrate 
that the hybrid method provides a much more useful solution than the perturbation solution 
alone. Some general comments about the method and a discussion of its possible application 
to certain problems associated with domain decomposition are discussed in Section 6. 
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2. Description of the Method 


The method we wish to describe is a two-step hybrid analysis technique based on the 
successive use of perturbation expansion methods and the classical Bubnov- Galerkin approx- 
imation technique. To illustrate the general ideas of the method, suppose we are seeking (an 
approximation to) the solution u to the problem 

(2.1) C(u,e) = 0, 

where C is some partial differential operator and e is a parameter. (Although we restrict 
our attention to two-dimensional elliptic boundary value problems in this paper, we believe 
that the method has a much wider range of applicability. Hence, we formulate the method 
in terms of a general operator £.) Here (2.1) holds in some domain V , and, in addition, 
must satisfy certain conditions on the boundary of V, which we denote by dV. Without 
loss of generality, we can assume that these boundary conditions are homogeneous in ti. 

In the first step of the method, we generate the coordinate functions in a perturbation 
expansion of u about one or more specific values of the parameter e, say about e - e p , p - 
1,2 , . . . , P. In the second step, we construct new approximate solutions consisting of sums of 
some of these perturbation coordinate functions, each multiplied by an unknown amplitude, 
and then determine these amplitudes by using the Bubnov-Galerkin method. 

To describe this idea in more detail, suppose that the solution to (2.1) can be expanded 

about e = £ p into a series of the form 

(2.2) u = U ^ P '^ a i( e ) + G(a5f ji+1 (e)), 

where {<*5(e)} is an appropriate asymptotic sequence of gauge functions and each u (pjl can be 
determined completely by standard perturbation methods (see, for example, Nayfeb [8]). By 
our assumption on the boundary conditions on u, each satisfies homogeneous boundary 

conditions on &D. 

A subset of all of the perturbation functions {u (p,i) } is now chosen as the set of coordinate 
functions for the Bubnov-Galerkin technique and approximations u for u are sought in the 

form 

N 

(2.3) «= 

j=l 

where the (unknown) parameters {^(e)} represent the amplitudes of the coordinate func- 
tions uW. Here each «« is one of the perturbation coordinate functions u^\ and hence 
ti satisfies the boundary conditions of the problem for any choice of amplitudes {S jiN }. To 
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determine these amplitudes, we apply the Bubnov-Galerkin technique to the governing equa- 
tion (2.1). Thus, we substitute (2.3) into (2.1) and require that the residual be orthogonal 
to the N coordinate functions over the domain V, i.e., 

( 2 - 4 ) / X £ “ (J) Sj ' N ( £ ^' e j uW dx = 0 > 1 <k<N. 

Equations (2.4) represent a set of N equations for the N unknown amplitudes. While (2.4) 
must, in general, be solved numerically, solving it is much simpler than numerically solving 
(2.1). In particular, for a fixed value of e, the solution to (2.4) is a point in W-dimensional 
space, where N is reasonably small, while the solution of (2.1) is a continuous function of 
the variables x. 

We should note that this particular choice of coordinate functions overcomes the main 
drawback of the Bubnov-Galerkin method, which is the difficulty, from a practical point 
of view, of selecting a small number of good coordinate functions. By the way they are 
constructed, the perturbation coordinate functions are (under certain assumptions) elements 
of a set of functions which span the space of solutions in a neighborhood of their point of 
generation. Thus, they characterize the solution u in that neighborhood. We also observe 
that, in many applications, the functions u( p,t ) are determined by solving a set of linear 
equations, even though the original operator C may be nonlinear. 

For the three applications discussed in this paper we use P = 1 only. 

3. A Simple Example 

To illustrate the basic features of the hybrid method, we consider the following simple 
two-dimensional example. We define u(x,y,e) as the solution to the problem 

(3-1) V 2 u + e sin(x) cos(y) u x + 2e sin(x) sin(y) = 0, (x,y)eV, 

(3-2) with u — 0 for (x,y)e3Z>. 

Here V = {(x,y) : 0 < x,y < tt} and dV denotes the boundary of V. In (3.1), the sym- 
bol V denotes the usual two-dimensional Laplacian operator, and the subscript x denotes 
differentiation with respect to x. The solution seems to be positive over the whole domain 
V for positive values of e. It is invariant under the 180° rotation x-> 7 r-x, y ->i r - y. 
Equation (3.1) is also invariant under the transformation e -v -e and u(x,y) -* -u(x,ir-y). 
Thus, we may focus our attention on solutions for positive e. Figure 1 shows surface and 
contour plots of the solution for three values of e. For large values of e, the solution exhibits 
a number of boundary layer properties (see Section 6). 
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Step One: For small values of e, we 


construct a regular perturbation expansion of u in 


the form 

( 3 . 3 ) u=-£e i uM(x,y) + 0(c N+1 ), 

j= 1 

where each of the perturbation coefficient functions is independent of e. To determine 
these functions, we substitute (3.3) into (3.1) and (3.2) and equate the coefficients of like 
powers of e. In this way, we are led to the following system of equations, from which the 
can be determined recursively. 


(3.4) 


V 2 u (,) 


J —2 sin(x) sin(y) if j — 1 

( — sin(x) cos(y) u^J~^ if 2 < j < N, 


( 3 . 5 ) u {j) = 0 for x = 0,7 r; and y = 0 ,tt for 1 < j < N. 

In particular, using (3.4) and (3.5) we find 

= sin(x) sin(y), ^ sin( 2 x) sin( 2 y), 


(3.6) 


u (3) _ jQj sin(3x) sin(3y) + Q) sin(3x) sin(y) 

- sin(x) sin(3y) - sin(x) sin(y)| , 

“ (4) = (49I52) s > n (^v) + (73^00) s * n(4x) s i n (2y) 

- (19555) sin(2x) sin(4y) - (^5) sin(2 I )»m(2 v ). 


It is easy to determine the general form of each u (fc) and to construct a recurrence relation 
for the numerical coefficients involved. Thus, many more terms in the expansion can be 
computed in a straightforward manner. 


Step Two: In the second step of the hybrid method, we use the perturbation coordinate 
functions {u^} obtained in Step One to construct new approximate solutions u in the form 

(3.7) u= X) 

n=l 

where N is a relatively small integer and the new amplitudes will be determined using 

the Galerkin technique. (We note that the u defined by (3.7) satisfy the boundary condition 
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(3.2) for any choice of the {<5^,#} because of the conditions (3.5) on the functions {u^}.) 
The procedure is as follows. We substitute (3.7) into (3.1) and require the residual to be 
orthogonal to each of the perturbation coordinate functions, i.e., 

J + € s * n ( x ) cos(y)u x + 2e sin(i) sin(y)] u (k) dx dy = 0, 

which may be written in the form 


(3.8) 

where 


N 


(3.9) 


"h ^ ^ij] ^ 1 ^ k < Ny 


c k,j = J J^y 2 u^]u^ dxdy, dk,j = J J^[sin(x) cos(y)u^]u^ dx dy, 

b, = -j j 2Mx) sin(y) dx dy . 


For this simple example, since the original problem is linear, the equations (3.8) to determine 
the {6j,w} are linear. Further, since the integrals in (3.9) can be evaluated exactly, it is 
possible in this case to evaluate the amplitudes as rational functions of c. However, we remark 
that for most applications (especially nonlinear ones) the evaluation of the amplitudes must 
be carried out numerically, typically by an iterative process for gradually increasing values 
of the parameter c. 

For N = 2, we use the expressions (3.6) to evaluate the coefficients {ckj}, {<4^} and 
{6j} explicitly and find 


(3.10) 


1 i .l/ioo) ^2,2 


1 + c 2 /128 : 


1 + e 2 /128’ 


In a similar manner, using the symbolic manipulation system Mathematica [15], for IV = 3 
we find 


(3.11) 

while for N = 4 


^ 1,3 — £> ^ 2,3 — 


1 + c 2 /60 


j ^ 3,3 = £ ^ 2 ,; 


, 3 > 


f _ 56524800c + 1216128 c 3 

1,4 “ 56524800 + 1216128 c 2 + 2141 c 4 ’ * 2 ’ 4 = e 5l>4 ’ 


(3.12) 


^3,4 = 


56524800 c 3 


56524800 + 1216128 e 2 + 2141 e 


7 > ^ 4,4 — £ ^ 3,4 
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Using Mathematica, we have obtained explicit expressions for the S jiN for N < 8. In each 
case, we find that 

(3.13) Sj.N = + 0(e N+1 ) as e -> 0, for 1 < j < N. 

Hence we have verified that our N-term hybrid solution (3.7) agrees with the N- term per- 
turbation solution (3.3) as e — * 0 through N = 8. 

In Figure 2, we have plotted the relative X 2 -error between numerical solutions to the 
problem (3.1 )— (3.2) and either hybrid solutions (3.7) or perturbation solutions (3.3) for dif- 
ferent values of N. Figure 2 demonstrates that: (1) there is a finite radius of convergence 
for the perturbation solutions; (2) the hybrid solutions are more accurate than the corre- 
sponding perturbation solutions, especially for larger values of e; and (3) for a fixed N the 
accuracy of the hybrid solutions gradually deteriorates as e increases. We discuss this exam- 
ple further in Section 6. The numerical solutions were obtained using a simple second order 
finite difference method (coupled with SOR) on a 41 by 41 uniform grid. For large values 
of £, the numerical solutions are more accurate than the perturbation and hybrid solutions, 
although for small e the reverse is true. 

4. Circular Cylinder in Slightly Compressible Flow 

We now apply our hybrid method to a problem involving a nonlinear partial differential 
equation. As a model of this type of problem, we consider the problem of determining the 
steady, two-dimensional flow of an inviscid, compressible, perfect gas past a circular cylinder 
of unit radius without circulation (see, e.g., Van Dyke [10]). If we introduce the usual polar 
coordinates (r, 6) and express the fluid velocity q in terms of a potential function <p as 

q = U v{ (r + r -1 ) cos(0) + <p(r, 9, M) j , 

we find that, for r > 1, ip satisfies 

C(ip, M) = V 2 tp-(±)M 2 {[(1 -r- 2 )cos(6) + cpr}Qr 

(4.1) 

-(-[— (r -1 + r~ 3 ) sin (8) + r~ 2 ipe\Qe + (7 — 1) Q = 

where 

Q = r -4 — 2 r~ 2 cos(20) -f 2 cos(0) (1 — r~ 2 )<p r — 2 sin(0) (r -1 4- r~ 3 )<ps 
+ <p 2 r +r- 2 <fl 

with w = 0 on r = 1, and <p = O (r' 1 ) as r -* oo. Here M = Ufc is the free stream Mach 
number, c is the speed of sound in the gas, and 7 is the adiabatic ratio. For small values 
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of M the flow streamlines are shown in Figure 3. For subsonic flow, the maximum velocity 
occurs at the surface of the cylinder at 9 = 7 t /2 and 0 = 37r/2. We now apply the hybrid 
method as outlined in Section 2 to this (nonlinear) problem, where M plays the role of the 
parameter £. 


Step One : For small values of AT, we use the regular perturbation method to obtain 
(4.2) (p(r t d,M) = Yl M2j ¥> (i) (r, 0) + 0{M 2N+2 ) as M - 0, 

J = 1 

where each satisfies a Poisson equation in the region exterior to the cylinder, satisfies the 
condition — 0 on r = 1, and is 0(r -1 ) as r — ► oo. The general form of the perturbation 
coordinate functions is given by 


(4.3) 


f {,) = fi.k{r) cos((2 k + 1)9], j > 1. 

k= 0 


In particular, we find: 

= (i) r ' 1 - (D r_1 + {ji) r ' 5 > 

< 4 -4) ^ = -(i) r * 1 + (n) r ‘ 3 - 

Van Dyke and Guttmann [12] have computed 29 terms in the expansion (4.2) using a FOR- 
TRAN program and have presented a rather detailed analysis of the convergence of the 
series as N — > oo. We have used the symbolic computation system Mathematica [15] to 
compute 5 of these terms in exact rational arithmetic, carrying 7 as a parameter, and 11 
terms for the special case when 7 = 7/5. Although the general form of each is known 
and only certain numerical coefficients need to be determined, the amount of computational 
“labor” (using either purely numerical or symbolic computation) increases significantly as 
N increases. Consequently, it is desirable to obtain as much information as possible about 
the solution from the first few perturbation coefficient functions { ip W}. This we do in Step 
Two below. 


Step Two : We now seek new approximate solutions (p in the form 

(4.5) <p( r , 0,M) = f>j, n(M) 0), 

j= 1 
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where the new “amplitudes” {$,-*} must be determined. To determine the {^}, we sub- 
stitute (4.5) into the governing differential equation (4.1) and require that the residual be 
orthogonal to each of the perturbation coordinate functions i.e., 


(4.6) 


ff £(<£, M) rdrd6 = 0, 


1 < Jfc < N. 


Equations (4.6) are a system of N (cubically nonlinear) equations to determine the JV am- 
plitudes {Sj t N(M)}. They must, in general, be solved numerically. However, this is straight- 
forward to do using a standard method such as Newton’s method. In particular, by starting 
at “small” values of M, where we expect S jt n « M 2j , and then proceeding to larger values 
of Af, the solution of (4.6) can be obtained in an efficient manner. 

For the special case when N = 1 , we let <p = S lt i<? (1) , where is given by (4.3) and 
(4.4) with j = 1, and we find that equation (4.6) yields the following cubic equation for 

(4 8 — M 2 [ 1 + ci 8 + C2 8 2 + C3 £ 3 ] = 0, 


where 

Cl = (4736 + 1111 7)/7532, c 2 = (147139 + 1582487)/1084608, 

(4.8) 

c 3 = (9182 + 815931 7)/23861376. 

From (4.7) and (4.8), we see that 

(4 9) 8 = M 2 + ci Af 4 + 0{M 6 ) as M -» 0, 

and hence our hybrid solution reproduces the first term of the perturbation solution as 
M -» 0. Also, the cubic equation (4.7) has one negative real root and two positive real roots 
for 0 < M < M c « 0.5389247. For M > M c , there are no positive real roots of this equation. 
In Section 6 we discuss a possible interpretation of M c in relation to the convergence of the 
series (4.2). 

In a similar manner, for N = 2, 3, 4, and 5, we have used Mathematica to perform the 
integrations appearing in (4.6) and have expressed the resulting equations in exact rational 
arithmetic. We then expressed the solutions {£j,iv} as Taylor series in M 2 and, in each case, 

we found that 


(4.10) 8 jiN = M 2j + 0(M 2N+2 ) as M — ► 0, fori <j<N. 

Thus, for these values of N we have verified that our hybrid solution (4.5) agrees with the 
first N terms in the perturbation expansion (4.2) as M -* 0. In addition, for each of these 
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values j>f N, we have also estimated the quantity M c , i.e., the value of M for which, if 

M > M c , Hie system (4.6) does not appear to have physically meaningful solutions. These 
values of M c are summarized in Table 1. 

As an application of the results we have obtained, we estimate the critical Mach number 
AT, associated with this flow. Here M* is defined as the value of the free-stream Mach number 

at which the flow locally becomes sonic for the first time and is determined as the solution 
of the equation 

( 4 - U ) M*=2/[(7 + 1) 9 Lx-7 + 1]. 

Here q max = 2 - d<p/dO evaluated at r = 1, and 9 = 7r/2 is the maximum surface speed. 
Van Dyke and Guttmann [12] have used 29 terms in the expansion (4.2), along with a 
battery of numerical techniques (see also Andersen and Geer [1] and Van Dyke [11]), to 
obtain the estimate M. = 0.39823780 ± 0.0000001. In contrast, using our one-term hybrid 
approximation we find that our approximation ^ for 9max is given by w = 2 + 
( 7 / 6 )^ i,i- Using equations (4.7) and (4.8) to compute for a particular value of M, we 
find from (4.11) that, for 7 = 7/5, our one-term hybrid solution yields the approximation 
M, = 0.407257966. This value differs from Van Dyke and Guttmann’s estimate by only 
about 2.3%. 

In a similar manner, we find approximations to both M c and M, using N = 2, 3,4, and 5 
terms in our hybrid approximation (4.5). The estimates we obtain, along with those obtained 
by using N terms in the perturbation expansion (4.2), are summarized in Table 1. 

TABLE 1 

Approximations to M* 


N 

Perturbation 

Hybrid 

M c I 

1 

0.420943 

0.407258 

0.5389 

2 

0.409239 

0.401704 

0.4990 

3 

0.404577 

0.399699 

0.4740 

4 

0.402274 

0.398948 

0.4606 

5 

0.400979 

0.398608 

0.4513 

6 

0.400187 



7 

0.399671 



8 

0.399319 



9 

0.399071 



10 

0.398891 



11 

0.398757 
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In Section 6, we discuss certain implications of the values which appear in Table 1. For 
the present, we just note that the hybrid approximation to M, based on N terms has about 
the same accuracy as the perturbation approximation based on 2 N terms. We have not tried 
to push our technique nearly far enough to begin to compete with Van Dyke’s results. We 
only observe that with a small number of terms we have a better approximation. 


5. Sink or Source Near a Free Surface 

We now wish to apply our hybrid method to a problem for which the nonlinearity enters 
through the boundary conditions; in particular, a problem in which the boundary itself is 
unkn own and must be determined as part of the solution to the problem. As a model of 
this type of problem, we consider the problem of determining the steady, two-dimensional 
shape of a free surface when a sink or source is placed above the surface. This problem has 
been studied using a variety of numerical techniques by several investigators, e.g., Tuck and 
Vanden-Broeck [9] and Vanden-Broeck and Keller [14], and also from a perturbation point 
of view by Vanden-Broeck, Schwartz and Tuck [13]. 

We consider a sink (or source) of strength Q located a distance h above the undisturbed 
height of a free surface. We introduce a rectangular coordinate system ( h x,hy ), with gravity 
acting in the negative y-direction and with the origin a distance h below the sink (see Figure 
4). We let the velocity potential be $ = Q { (^;) logt® 2 + {y ~ l) 2 ] + <p( x > 3/)} an( l denote the 
elevation of the free surface by y = tj(x). Then the problem to determine <p and rj becomes 


(5.1) 


VV = o, y > v(. x )> 


with 

(5.2) 


Bifay) = v’y - VxV - 



_ Q 

x 2 + fo - l) 2 


and 

(5.3) 


B 2 (v,ri) = V ~ + V 2 ) 5 -£(1 + t) 2 ) 


<Px + 


( 2 ^) x 2 + ( 7 ? — l) 2 


= 0 


on y = j?(x). Here = d/dx, £ = Q 2 /2gh 3 , and a = 2^h/pQ 2 , where g is the acceleration 
due to gravity, p is the density of the fluid, and qr is the surface tension coefficient. Equation 
(5.2) is just the kinematic boundary condition on the free surface, while equation (5.3) follows 
from Bernoulli’s equation and the condition on the jump in pressure at the free surface due to 
surface tension. We are interested in finding approximations to the solutions <p = <p(x,y, e) 
and 77 = r)(x,e) to equations (5.1)— (5.3) for small values of e. 
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Step One : For small values of e, we look for regular perturbation solutions for and 77 
in the form 


(5.4) 


'f = Z) t3 V> (j) + 0{e N ), V = J2 £j Vj + 0(e N+1 ) as e -* 0. 
3=0 j = 1 


Using the expressions (5.4) in equations (5.1)-(5.3), we find that the terms which are 0(1) 
in e in equations (5.1) and (5.2) yield the result 

( 5 -5) ipW = (-^- ) log(x 2 + (y + l) 2 ). 

Then the terms which are 0(e) in equation (5.3) yield the expression 

( 5 - 6 ) Vi = (^j) a: 2 (l + x 2 ) -2 . 

In a similar manner, the terms which are 0(e) in equations (5.1) and (5.2) yield the result 


(5.7) 


<pM = 


3(x 2 — (y + l) 2 ) , 2(y + l)[(y + l) 2 — 3 x 2 ] 
8 tt 3 [ (x 2 + (y + l) 3 ) 3 + (x 2 + (y + l) 2 ) 3 

and the terms which are 0(e 2 ) in (5.3) yield the expression 


(5.8) 


Vi = 


3 

x 2 (l — 6 x 2 + x 4 ) 

2a 
1 

1 — 8 x 2 + 3 x 4 

2tt 4 

(1+X 2 )5 

1 2 
7T^ 

(1 + X 2 ) 4 


Continuing this procedure with the aid of Mathematica, we find that the general form of 
the perturbation coefficient functions </jU) and rjj , j > 1, is 


(5.9) 


3 j 


<p M (x,y) = (l/n 2j+1 )j^c jtk iJ>W(x,y), 


k=2 


J-l 


Vj(x) = ak PjA x 2 ) (i + x 2 ) (3i 1_A:) . 

k=0 

Here each c hk is a constant, each pj <k is a polynomial in x 2 , and the functions are 

defined by 

v ,(0) = log[x 2 + (y + 1) 2 ]5, A k) = (d/dy)^- 1 ) for k > 1. 

In particular, we find 

(x) _ x 2 (45 — 1909 x 2 + 6890 x 4 — 2786 x 6 + 9 x 8 — 9 x 10 ) 

V3[X 32 7T 6 (1 + x 2 ) 8 


a (12 - 853 x 2 + 4008 x 4 - 2658 x 6 + 148 x 8 - x 10 ) 

4 7T 4 ( 1 + X 2 ) 7 

24 a 2 (2 — 33 x 2 + 40 x 4 — 5 x 6 ) 

7r 2 (l+x 2 )6 
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C 2>2 = -c 2>3 = (9 — 87r 2 a)/128, c 3>4 = (27 - 807r 2 a)/384, 

C2 5 = (3 - 47r 2 a)/64, c 2 . 6 = 3/640. 

From the analysis of this perturbation series with zero surface tension by Vanden-Broeck 
et al [13] (see Section 6 for more discussion of their results), our formal senes expansion of 
the solution for this example does not converge for any non-zero value of e. However, espi e 
this lack of convergence, we proceed to Step Two of our hybrid method. 

Step Two : Using the expressions above for the perturbation coordinate functions v) 

and r 7 j(x), we look for new approximate solutions in the form: 

(5.10) f 7 (x,e) = > £(*.»»«) = 

3 = 1 J 

where the new “amplitudes” {«, „} must be determined. (We note that <e satisfies V\> = 0 
for any choice of the amplitudes {«*»}.) To determine these amplitudes, we substitute the 
expressions (5.10) into the boundary conditions (5.2) and (5.3) and require that the residuals 
be orthogonal to appropriate sets of test functions {a*} and {ft}, 

(5.11) /”Bi(v, «)<**(*) *>= 0. £>(lM) «*)■<* = »■ ^ k ± N - 

Once the test tactions have been selected, equations (5.11) are a set of 2 N equations 

for the 2N unknowns {«,.»,, = 1,2 2W}. For simplicity we select as test functions 

ai = ft = tji(i). The resulting hybrid solutions are discussed in the next section. 

6. Discussion 

The hybrid perturbation-Galerkin method, as we have described it here, is a semi- 
analytical, semi-numerical technique which appears to have the potential of being a useful 
tool both to complement and to supplement existing standard methods of analysis is 
semi-analytical in the sense that some of the analytical structure of the solution is deter- 
mined by first constructing one or more perturbation approximations to the solution is 
semi-numerical in that new amplitudes of the perturbation coordinate functions typica y are 
determined numerically by solving standard Galerkin equations associated with the prob- 
lem The method seems to have considerable potential to be an effective problem-solving 
tool because the perturbation coordinate functions appear to be very effective trial tactions 
in the Galerkin approximation. Intuitively this is because, by the way they are constructe 
they describe the basic nature of the solution, at least for parameter values near their point 

of generation. 
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In our first example (Section 3), we analyzed a simple, lineax boundary- value problem 
on a bounded domain, for which an arbitrary number of terms in the regular perturbation 
series can be computed in a straightforward manner. By examining the singularities of the 
amplitudes in our hybrid solution, it appears that the convergence of the perturbation 

series is limited by a pair of complex conjugate singularities located near or at ±6.70z. (This 
follows by analogy from the analysis of several similar problems, involving ordinary differ- 
ential equations, for which the singularities of the {6 jt in the complex e-plane converged 
to the singularities of the actual solution as the number of terms in the approximation in- 


creased. See especially [2]. In fact the location of the poles of the { S jiN } which are closest 
to the origin form the sequence ±7.75 i, ±7.15 i, ±6.803 i, ±6.721 i, ±6.7048 i, ±6.70283z, .... 
Another pair of singularities is indicated near e « ±9i .) These singularities appear to have 

no immediate interpretation, but nonetheless limit the direct usefulness of the perturbation 
expansion. 

In Figure 2, we have plotted the relative Z 2 -error between numerical solutions to problem 
(3.1)-(3.2) and perturbation or hybrid approximations. The figure clearly illustrates that 
the hybrid method allows us to obtain useful information about the solution well beyond 
the radius of convergence of the perturbation solution. Also, there is no indication that the 
range of convergence of the hybrid solution is limited. 


In Figure 6, we have plotted the values of the solution at the center value x = y = 7r/2 as 
calculated numerically and by the use of the hybrid technique, where the number of terms 
(N) ranges from 1 to 8. It is interesting to note that our hybrid solutions appear to form an 
alternating sequence which brackets the true center point solution. Figure 1 shows that the 
solutions have their peak ^elevations” at the center point. 

The contour plots shown in Figure 7 demonstrate that higher-order perturbation func- 
tions can be very similar to one another. This suggests that the perturbation functions which 
are of higher order than these shown may be fairly well approximated by linear combinations 
of the functions shown (see observation (2) in the introduction). Figure 7 also demonstrates 
that the perturbation functions may have higher symmetries (group order = 4) than the so- 
lutions (group order = 2) which they are used to approximate. Such additional symmetries 
may be exploited in the steps of the hybrid method which are computed numerically. The 

symmetries of the perturbation functions follow a regular pattern which can be determined 
a priori. 

From a singular perturbation point of view, as e becomes large the solution to problem 
(3.1)-(3.2) takes on a somewhat complicated form. The leading term in the outer expansion 
of u is 2 (t r - x)tan(y) for 0 < x < tt, 0 < y < tt/ 2, and is -2*tan(y) for 0 < x < 
tt, tt/2 <y < tt. Boundary layers develop along x = 0 for 0 < y < tt/ 2 and along x = tt 
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for 7r/2 < y < tv. In addition, an internal layer forms along y = 7r/2, as well as additional 
boundary layers at (0,7r/2) and {tv ,tv j 2). From the contour plots shown in Figure 1, it may 
be seen that the boundary and internal layers are beginning to be formed by the hybrid 
approximations as e increases, even though these approximations are based on the small c 
expansion of the solution, where no such layers are evident. 

This last observation suggests that the hybrid method might be a useful tool for the 
general problem of determining an appropriate decomposition of the domain for a purely 
numerical solution to a specific problem. In particular, suppose that some “singular” behav- 
ior of the solution is anticipated for large values of a parameter appearing in the problem 
formulation. To help determine the location and, to some extent, the nature of this singular 
behavior, the hybrid method could be applied in the following way. First, a (regular) pertur- 
bation expansion of the solution is constructed for small values of the parameter. The hybrid 
method is then applied, using the small parameter perturbation coordinate functions, and 
then the behavior of the hybrid approximation is noted as the parameter value is increased. 
From the example of Section 3, as well as several other examples we have examined in some 
detail, it appears that the hybrid solution simulates at least some of the singular behavior 
of the exact solution as the parameter value is increased, and that this simulation becomes 
more accurate as the number of terms in the approximation is increased. 

In our second example (Section 4), which involves an infinite domain, the regular per- 
turbation expansion of vp is straightforward, although tedious, to compute. Frankl and 
Keldysh [4] have proved that <p is analytic in M 2 about M 2 = 0 and hence the pertur- 
bation expansion (4.2) is actually the Taylor series expansion of the exact solution. They 
did not, however, determine the radius of convergence M c of the series. Through a care- 
ful analysis of the first 24 terms in the series, Van Dyke and Guttmann [12] report that 
M c = 0.402667605 ± 0.000000005, which is about 1.11% greater than their estimate of M*. 
If this were indeed the case, it would imply that shock free solutions exist for a continu- 
ous range of values of M above the critical Mach number. However, Van Dyke [private 
communication] now believes this result to be in error and that M c = M„. 

Although Van Dyke and Guttmann were unable to determine the nature of the singularity 
which limits the convergence of the series, we can safely assume that the singularity is related 
to the formation of a shock and hence represents a reed singularity with a definite physical 
interpretation. The numbers M c presented in Table I are estimates for the value of M above 
which no physically meaningful hybrid solution to the problem exists. In this sense, they can 
be thought of as estimates of the radius of convergence of the series (4.2). If more of these 
estimates were computed and if it could be shown that they converge to the same limit as 
the estimates M ,, this would lend support to the conjecture that M c = M t . 
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The solution to our third example (Section 5) involves two unknown functions (i.e., 
the potential function and the shape of the free surface), and the associated perturbation 
expansions are somewhat more difficult to compute than the expansions associated with 
the first two examples. However, the general form of the terms in each expansion can be 
determined and, in principle, only certain numerical constants need to be computed. Using 
a complex variable formulation of this problem (with zero surface tension) Vanden-Broeck et 
al [13] numerically computed these constants in the first 34 terms in each expansion. They 
analyzed the resulting series for the free surface and found it to be of exponential-integral 
character (i.e., the coefficient of e n behaves like n! for large n). Consequently, the series 
diverges for all nonzero values of e. They use a number of techniques to “sum” this divergent 
series and show that this analysis leads to a free surface shape with a jump discontinuity 
due to the location of a branch cut. One such free surface profile corresponding to our 
£ = 7r 2 /5 = 1.974 is indicated by a dashed line in Figure 5. Our hybrid solution with N = 2 
is shown by a solid line. Vanden-Broeck et al indicate that the jump discontinuity in their 
solution could be removed a applying a certain iterative procedure, but the corresponding 
profile was not displayed. Using a numerical, series truncation method on another complex 
variable formulation of the problem, Tuck and Vanden-Broeck [9] found a cusp like solution, 
corresponding in our notation to a value of e = e c « 6.311. Vanden-Broeck and Keller 
[14], using a similar method, treated the same problem we are considering, except that a 
horizontal fixed surface is present at a finite distance above the sink. They found solutions 
corresponding to our problem for values of e larger than e c , but did not find any steady 
solutions for 0 < £ < £ c . They indicate, however, that solutions with waves may exist for £ 
in the range 0 < £ < e c . 

The hybrid approximations presented in Section 5 appear to give “physically realistic” 
approximations to a steady solution to the problem, as indicated in Figures 4 and 5. In 
[2] we applied our hybrid method to an eigenvalue problem associated with an anharmonic 
oscillator. For this problem, the perturbation series also has a zero radius of convergence, 
but the hybrid approximations converge monotonically to the true solution as N increases. 
Based upon this experience, it is certainly possible that our hybrid approximations to the 
free surface for values of £ in the range 0 < £ < £ c are converging to a steady state solution. 
Obviously, this difficult problem needs much more investigation, e.g., either the development 
of a numerical method which will yield the steady solution for £ in the range 0 < £ < £ c , or 
a proof that no such solution exists in this range of parameter values. 
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FIGURE CAPTIONS 


Figure 1. Surface and contour plots for 8-term hybrid and numerical solutions of the simple 
partial differential equation of Equations (3. 1)— (3.2). The hybrid solutions shown for £ = 1 
and £ = 5 are highly accurate. Some irregularities may be seen in the relatively flat portions 
of the hybrid solution for e = 20 which do not appear in the numerical solution for £ = 20. 
Nevertheless, the relative .La-error is less than 1%. 

Figure 2. Relative La-errors of perturbation and hybrid solutions of the simple PDE problem 
as a function of the parameter e. The number of terms N ranges from 1 to 8. The log-log 
plots of the perturbation errors fall nearly on straight lines which intersect one another near 
the radius of convergence, R = 6.7026. The hybrid errors for given N and £ are significantly 
less than the perturbation errors for the same N and e and show no signs of divergence. 

Figure 3. Streamlines for two-dimensional flow around a cylinder for Mach number M = 0, 
the low velocity limit of the solutions to Equation (4.1). 

Figure 4. Free surface profile for the two-dimensional flow induced by a sink (black dot) 
above a liquid surface (see Equations (5.1)-(5.3)). The effects of surface tension are ignored. 
The horizontal and vertical scales differ. 

Figure 5. Free surface profiles for £ = 7r 2 /5 with zero surface tension as computed by the 
hybrid method (solid line) and as computed by Vanden-Broeck et al (dashed line). 

Figure 6. Maximum values of the solution (x = y = 7r/2) of the simple PDE problem 
of Section 3 as computed numerically (dots) and by the hybrid method (solid lines) using 
from N = 1 to N = 8 terms. The results for N = 8 are useful well beyond the radius of 
convergence 3 R } of the series expansion. 

Figure 7. Contour plots of the first eight perturbation functions for the simple PDE problem 
of Section 3 demonstrate that the higher order perturbation functions can be very similar to 
one another and can have higher symmetries than the solutions they are used to approximate. 
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